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Abstract 

In this paper we discuss a new approach to the quasi-normal mode problem in general 
relativity. By combining a characteristic formulation of the perturbation equations with the 
integration of a suitable phase-function for a complex valued radial coordinate, we reformulate 
the standard outgoing-wave boundary condition as a zero Dirichlet condition. This has a 
number of important advantages over previous strategies. The characteristic formulation 
permits coordinate compactification, which means that we can impose the boundary condition 
at future null infinity. The phase function avoids oscillatory behaviour in the solution, and 
the use of a complex radial variable allows a clean distinction between out- and ingoing waves. 
We demonstrate that the method is easy to implement, and that it leads to high precision 
numerical results. Finally, we argue that the method should generalise to the important 
problem of rapidly rotating neutron star spacetimes. 



1 Introduction 

The dynamical oscillations of compact objects is a problem of great relevance for general relativistic as- 
trophysics. Interesting questions range from observational to theoretical, with the possibility of detecting 
gravitational waves from pulsating neutron stars and black holes and using the signals to infer source param- 
eters providing strong motivation for detailed studies. From the more theoretical point of view, the stability 
properties of these astrophysical bodies are intimately linked to their oscillation spectra. For rotating neu- 
tron stars, this is clearly demonstrated by the gravitational-wave driven instability that was first discussed 
by Chandrasekhar ,1^, Friedman and Schutz [2]. 

In general relativity, non-radial fluid oscillations radiate gravitational waves. Hence the oscillation mode 
problem is conceptually different from that in Newtonian gravity. In order to determine the pulsation modes 
one must impose an outgoing-wave boundary condition at infinity. It is well documented that this leads 
to technical difficulties if the quasi-normal modes are rapidly damped. The main reason for this is that on 
a spacelike hypersurface, as obtained by assuming that the perturbation quantities have a exp{iLut) time- 
dependence, the required solutions grow exponentially towards infinity. In order to impose an accurate 
outgoing- wave condition one must be able to filter out the exponentially decaying ingoing- wave component. 

Various methods have been devised to handle this problem, see the review by Kokkotas and Schmidt [5] 
for an exhaustive discussion. One approach, that has the advantage that it is relatively easy to implement, 
is based on the use of analytic continuation and integration of the perturbation equations for complex values 
of the radial variable. The method, which was first used to calculate accurate quasi-normal modes of black 
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holes 1] , has been successfully applied to the stellar oscillation problem [S] . 

The main outstanding problem in this area concerns the oscillations of rapidly spinning relativistic stars. 
In this case the perturbation equations that need to be solved in the exterior vacuum are no longer describable 
by a single wave-type equation (essentially because the spacetime is no longer of Petrov type D [5). This 
makes the solution more involved and, in particular, any method that relies on the solution of a single 
separated differential equation no longer applies. This problem has yet to be overcome. The best results 
that we have correspond on the one hand to the so-called neutral modes [71 [5] , identifying the points where 
the fundamental f-mode becomes susceptible to the Chandrasekhar-Friedman-Schutz instability, and on the 
other hand to modes determined after ignoring the metric perturbations (within the relativistic Cowling 
approximation) [9]. In the latter case one can use multipole formulae to estimate the damping/growth rate 
of the modes. Nevertheless, the situation is not satisfactory. After all, it is not clear how accurate the 
Cowling approximation will be for the various classes of modes that one may be interested in. There is a 
clear need for a more detailed solution to the mode-problem for rapidly rotating neutron stars. 

This paper introduces a new strategy to the problem. The main idea is to use a characteristic formulation 
of the perturbation problem [10] to avoid the diverging eigenfunctions that plague the standard analysis. By 
considering two simple model problems, we will demonstrate the promise of this approach. Most importantly, 
we will show how one can formulate the outgoing-wave condition as a Dirichlet condition for a suitably 
introduced phase- function [llj . Taken together with the fact that a characteristic formulation permits 
coordinate compactification without the loss of resolution [12] , this enables us to impose the desired boundary 
condition at future null infinity with extreme precision. We illustrate the key features of the new method by 
applying it to the well-studied problem of axial spacetime modes of a uniform density star. Using a spectral 
approach to solve for the relevant phase-function in the exterior of the star, we obtain results with very 
high numerical precision. The reliability of the method is shown by considering an ultra-compact star, with 
R/M — 2.26, which has both very slowly and very rapidly damped modes. We confirm previous results for 
the slowly damped modes, and even provide some new results by identifying several previously unknown 
interface w-modes. This result is likely irrelevant from the astrophysics point of view, but it is conceptually 
interesting since it hints at the possibility that there may exist an infinite number of such modes. 

The numerical results that we report are in themselves not particularly interesting. The main achievement 
is the successful new formulation of the quasi-normal mode problem that should be straightforward to extend 
to the rapid rotation case. The characteristic formulation of the Einstein equations is of course well-known, 
and we cannot see any real difficulties in generalising our method to the resultant coupled perturbation 
equations. In fact, the possibility of having a Dirichlet condition at infinity in the computationally more 
involved two-dimensional problem is very attractive. Some details obviously remain to be worked out. The 
main issue concerns phase-functions for coupled equations, a problem that has already been considered in 
quantum scattering problems 13]. Hence, we feel very optimistic. For the first time we have a truly promising 
strategy for dealing with one of the most challenging problems in this area of research, the determination of 
oscillation modes and the associated gravitational- wave damping/growth rate for fast spinning stars. 



1.1 Problem setup 

As a model problem we will consider the mode problem in static spherically symmetric spacetimes. In 
particular, the wave equations we will discuss in this paper can all be written on the form 

{VaV -U)^P = (1) 

where Va is the covariant derivative on the submanifold orthogonal to the spherical symmetry surfaces (a 
runs from to 1) and U is an effective potential depending only on the radial coordinate. For vacuum 
perturbations the potential is of the form 

U='l±»-'-Eu, (2) 
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where 



U- = 1 (3) 
_ [l\l + lY-Ay + l2M{r-M) 
+ [{l-\){l + 2)r + QMY ^ ' 

for axial and polar perturbations, respectively. Here M is the mass, r is the Schwarzschild radial coordinate 
and Z > 2 is the usual angular "quantum number" originating from the separation of the angular dependence. 
Asymptotically the polar piece approaches a constant; 

l{l + l) + 2 
l{l + l)~2 

and it is easy to show that | < U+ < 2 everywhere outside r = 2M. It is evident that the asymptotic 
behaviour of the axial and polar perturbations is very similar. For this reason we restrict our attention, from 
now on, to the axial sector. 

We shall later, as a test case, consider axial perturbations of perfect fluid stars. These are also governed 
by an equations of the form ([T]) , with U given by 

1(1 + 1) 6m K . . 
U=^^~ — + ^ip-p) (6) 

where m — m{r) is the mass inside r, p is the energy density and p is the pressure [141 115] . 

In our analysis of the wave equation ^ we shall use our freedom of choice of coordinates. It is conventional 
to adapt the coordinates to the timelike Killing vector on the background and choose the remaining radial 
coordinate to be either the Schwarzschild area radial function r or the Regge- Wheeler tortoise radius . 
The former has the advantage of being invariantly defined as giving the area of a spherical symmetry surface 
through Area = 47rr^. On the other hand, the tortoise coordinate is intimately related to the characteristics 
of the wave equation for massless fields and is related to r (in vacuum) through 

n = r + 2M ln[C(r - 2M)] (7) 

where C is a constant related to translations. Since we are primarily interested in the asymptotics (which 
does not depend strongly on whether r or r, is used) we shall use r in the following. 

An alternative to using the Killing time as a coordinate is to use characteristic coordinates as introduced 
by Bondi et al.fW and Sachs [I7|. The relation between the different coordinates is 

x = r (8) 
u^t-r^r) (9) 

where u is the retarded time and the characteristic radial variable x is actually the same as the area radius 
and we shall soon denote it by r. However it is important to remember that 



d dx d du d d dr* d 
dr dr dx dr du dx dr du 

For vacuum spacetimes and for the coordinates discussed above the line elements can be written 
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(10) 



where 



ds^ = dsj +r^(dr +sin^6'd0^) (11) 

As] = -e^^du^ - 2dudr Bondi-Sachs (B) (12) 
ds] = -e^^dt^ + e-^^dr^ Schwarzschild (S) (13) 
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We have set a; = r and introduced the notation 



e^" = 1 14 

r 

The wave equation ([T|) can now be cxphcitly written down 

-2V',„r + e^-'V.rr + — V'.r " C/?/' = (B) (15) 

-e-'^V,** + e^^'i^.rr + '^i^.r - ^ (S) (16) 

Note that, in the characteristic formulation, there is no second derivative with respect to the "time" u. To 
find mode solutions we take the time dependence to be given by "0 oc e'""^ where T = t,u depending on 
coordinates. We let primes denote derivatives with respect to the spatial coordinate (remember that the 
primes are different due to (fTU]) !). 

e2>" + 2 - tio^ 'ip' - C/V^ = (B) (17) 
+ 7^ V'' - (f/ - e-^^'w^)^ ^ (g) (^8) 

We may note here that the algebraic transformation tp e^^'^^'ip brings the Schwarzschild equation to the 
form of the Bondi-Sachs equation. Thus for the simple case of static spacetimes considered here we can 
mimic the characteristic approach by this simple transformation. However, the case we are really interested 
in is the two-dimensional rapidly rotating problem where such a trick is not likely to exist. 



2 A pedagogical toy problem 

Since the main interest in this paper is the behaviour of wave solutions in the asymptotically flat region of 
spacetime far away from the source we will begin by considering the very illustrative toy problem in which 
M is set to zero. It should be intuitively clear, and it will later be demonstrated, that this problem has the 
same properties near infinity as the M 7^ case. A great advantage is that it allows for an exact solution. 
In fact, in Schwarzschild coordinates 

V> = V^[a„i?/+l/2(^r) + CoutHll\^^iLur)] (S) (19) 

where i?^^-* and iJ*^^^ represent the in- and outgoing Hankel (or Bessel of the third kind) functions and the 
constants refer to the proportion of out- and ingoing waves as can be seen from the asymptotic behaviour 

V'e'"* ~ Cine''^(*+'') + Coute*"^*-''^ (S) (20) 

As noted previously, in Bondi-Sachs coordinates the solution is just given by multiplying the Schwarzschild 
solution by e*"'"* = e''^''. In order to illustrate the nature of these solutions we will now examine the special 
case of Z = 2. Then the solution becomes, after some redefinitions of the constants, 

uj^r^ + Siujr ~ 3 cj^r^ - 3iu;r - 3 
V = Cin 5-^ e -Gout 5-^ e (S) (21) 

For definiteness we will plot the solutions for an arbitrarily chosen frequency co = 1 + i/10 in units such that 
the inner boundary is located at r = 1. The ingoing and outgoing parts of the solution in Schwarzschild 
coordinates arc plotted in figures [1^ and[T]D, respectively. If we imagine finding these solutions numerically 
we immediately note two problems. First the solutions are oscillating. This is clearly a problem since 
we are really interested in imposing boundary conditions at infinity and hence would like to compactify 
spacetime in order to cover it on a finite grid. It is evident, regardless of the number of grid points, that 



4 




10 15 20 25 30 





15 20 



Figure 1: The left set of panels (a-d) display the real (solid) and imaginary (dashed) parts of the wave 
function V' as a function of r. The right set (e-h) shows the phase function g{r). In each set of four panels 
the top two (a-b and e-f) show the solution in Schwarzschild coordinates whereas the bottom two (c-d, g-h) 
display the characteristic solutions. The left panels (a,c and e,g) show the ingoing solutions while the right 
panels (b,d and f,h) instead display the outgoing solution. In panels f and h we additionally plot a generic 
solution with Cm = 0.57, Cout = 1-3. See the main text for further discussion. 



the oscillating behaviour will not be resolved at some finite radius. The other problem is that the ingoing 
solution is exponentially decreasing outwards. Hence, it will be numerically impossible to distinguish the 
outgoing solution from a solution contaminated by an ingoing piece. In figures [TJ; and [TJi we instead plot 
the characteristic/Bondi-Sachs form. The outgoing solution is clearly much better behaved here and would 
in principle be possible to track numerically on a compactified grid. However, the problem with the rapid 
decay of the ingoing solution remains. 

The oscillation problem is often resolved by introducing a phase function. Several options exist which 
are essentially equivalent, see e.g. Here we simply define the phase function to be 

9 = T 22 

■0 

as suggested from WKB-type arguments. In figure[T]3-h we show the behaviour of the phase function. Panels 
[T^ and [if show the ingoing and outgoing solutions in Schwarzschild coordinates whereas [T^ and [iji arc the 
Bondi-Sachs forms. In panels [if and[Tfi we additionally show a generic solution with coefficients arbitrarily 
chosen to be Cin = 0.57 and Cout = 1-3. We see that although the phase function solves the problem of 
oscillating solutions we still cannot distinguish the outgoing solution from a generic one at infinity. The 
reason is that the transformation to the variable g is non-linear and pieces coming from the ingoing solution 
will always be exponentially suppressed for positive imaginary parts of w. It is clear that we must find a 
way to remove the exponential decay of the ingoing solution in order to find a successful numerical scheme. 
Before turning to this problem it is useful to note that although all solutions are asymptotically constant, 
the outgoing Bondi-Sachs version has the advantage of being asymptotic to zero. Hence, in that formulation, 
the outgoing wave boundary condition can be replaced by a much simpler zero Dirichlet condition. 
We turn now to the decay of the ingoing solution. The problem comes from the asymptotic factor 




Figure 2: Same as figure [T] but with g as radial coordinate. Hence, these solutions are obtained from 
integration in the complex r-plane. 

There is a neat trick to resolve this problem which has been used for quasi-normal modes, see e.g. [HIS]. The 
idea is to analytically continue the support of the dependent variable g or ijj and let the radial coordinate 
take on complex values. In this spirit we write 

r = B.+ ge'"^ (24) 

where R is the (real) radius outside which we arc interested in the solution and g and § are real. We now 
wish to choose the angle d in the complex plane such that the decay of the ingoing solution is avoided. 
Substituting ([^^ into we see that we must choose 

?9 = - arg(w) (25) 

Having made this choic^ we plot the various solutions again in figure [2l Note that all divergent/decaying 
behaviour is gone. We see that, in principle, there exist three possibilities to implement a successful numerical 
scheme. The characteristic approach aided by the complex radial coordinate as displayed in panels[2t-d allows 
distinction of the in- and outgoing solutions together with compactification and imposition of boundary 
conditions at or near infinity. The same applies to the phase function approaches as shown in figures [2l3-f 
(Schwarzschild) and[2fe-h (Bondi-Sachs). However, the Bondi-Sachs phase function has the advantage of going 
asymptotically to zero. Hence, in this approach, the outgoing-wave boundary condition may be replaced by 
a zero Dirichlet condition on the phase function at null infinity. This is of course a great simplification when 
extensions to more general circumstances (such as rapidly rotating sources) are considered. 

2.1 A numerical example 

In order to show that the approach can be accurately implemented we now solve a toy quasi-normal mode 
problem. As above we take M — and choose units in which i? = 1. To complete the problem we choose 
boundary conditions such that, in Bondi-Sachs coordinates 

^'\r=R = g\r=R^ = (26) 

^Obviously, a; = is a degenerate case in this respect and cannot be treated with the present approach. However, since the 
oj = problem for rotating stars has aheady been studied in detail by Stergioulas and coworkers [7j [S] this is not a serious 
drawback of our method. 
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I 


Frequencies, exact 


numerical 


2 


2i 


0.000000000 + 2.000000000i 


3 


i(\/5 + 5i) 


1.118033989 + 2.500000000i 


7 


5.071731535 + 3.511061460i 


5.071731535 + 3.511061460i 




2.947010505 + 4.736313493i 


2.947 + 4.74i 




0.9710593605 + 5.252625048i 


0.97106 + 5.2526i 



Table 1: Quasi-normal modes of the toy problem for various I. All the digits shown have converged and 
agree with the exact solution. In some cases we reach machine precision, see table [2] Note, however, the 
poor convergence for some of the I = 7 modes. See the main text for an explanation of this. 



whereas, as shown above, g satisfies a zero Dirichlet condition at infinity. The inner boundary condition 
is chosen purely for mathematical convenience. In order to be able to enforce the boundary conditions at 
infinity we compactify spacetime by introducing the radial coordinate x through 

X - ^ (27) 

ri + g 

where the real constant ri makes it easier to keep track of the dimensions. Note that 

r = R ^ g = Q ^ x = l (28) 
r = oo ^ g = oo ^ x = —I (29) 

so that the exterior spacetime is mapped to the domain — 1 < a; < 1. The reason we chose this range is that 
we base our numerical code on a Chebychev pseudo-spectral method, e.g. 18J. Using equations (|15p and 
([22|) we find that the phase function has to satisfy 

= 2 _ 4za;r2 _ 27-2/(^ + 1) 

(x + l)^^ (x + l)2^ [i?(l + a;)+r2(l-x)]2 ^ ' 

where r2 = rie'''. As will be seen in the next section, the irregular singular point at infinity {x = — 1) 
does not pose a problem for the outgoing solution since g ^ {1 + x)'^ in the neighbourhood of null infinity. 
The purely ingoing solution is (of course) also well behaved since it asymptotically cancels the two first 
terms on the right hand side of ()30|) . When we enforce the boundary conditions at infinity we add the point 
x = —1 and put g — exactly there. This is one of the advantages of the pseudo-spectral approach - we are 
actually able to enforce the boundary conditions at infinity. The family of solutions to ([3D]) then provide a 
function gs{^) — g{x = l,w) whose zeros correspond to the quasi-normal modes. At a practical level, we 
solve equation fSn]) by a Chebychev pseudo-spectral code using a Newton-Kantorovich scheme to turn the 
problem from a non-linear equation to iteratively solved linear equations (see e.g. the book by Boyd |18j for 
details on spectral methods). We iterate until the corrections to the solution are everywhere less than 10~^^. 
We solve the equation for a large set of frequencies. A contour plot of the absolute value of gs{'-^) is then 
used to locate the points in the complex frequency plane where 5s (w) is close to zero. These points are then 
used as initial guesses in a Miiller type root finder. We decided to use a spectral code since we expect this 
approach to generalise to the coupled partial differential equations which need to be solved in the rapidly 
rotating case. 

We examine the cases / € {2,3,7}. The eigenmodes in the positive quadrant are given in table [T] For 
the case I = 2 we see that we have the extreme case of a "mode" with purely imaginary frequency. In 
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N 


Frequency 


2- 


u;\ 


iTT - arg( 




40 


5.444541 X 10 


-3 _ 


f 1.997613i 


-2.379953 x 10 


-3 


-2.725517 X 10 


~3 


80 


4.063484 X 10 


-7 


f 2.000000i 


4.089420 X 10 


-8 


-2.031742 X 10 


-7 


120 


1.893739 X 10- 


11 


f 2.000000i 


-1.279776 X 10- 


11 


-9.468648 x IQ- 


12 


160 


3.011046 X 10- 


15 _ 


f 2.000000i 


1.332268 X 10- 


15 


-1.554312 X 10- 


15 


200 


-4.551599 X 10" 


16 , 


f 2.000000i 


4.440892 X 10- 


16 


2.220446 X 10" 


16 


240 


9.935105 X 10- 


16 - 


f 2.000000i 


-6.661338 X 10- 


16 


-4.440892 X IQ- 


16 



Table 2: Convergence test for the I = 2 quasi-normal mode of the toy problem. The two right-hand columns 
show the error in the numerically obtained absolute value and phase of the frequency as compared with the 
exact value to = 2i. Machine precision is reached exponentially fast as the resolution (N) is increased. 



table [2] we show the convergence of the solution as a function of grid points N. As we can see machine 
precision (in this case double ^ 10~^^) is reached exponentially fast, a trade mark of a spectral code. The 
convergence is equally impressive for I ~ 3. For the case I = 7, however, we run into problems. The first 
mode K, 5.07 -|- 3.51z) is found without problems, but as the argument of the frequency starts growing the 
code runs into convergence problems. The reason for this behaviour can be seen from the exact solution. 
For any I > 2 the solutions have the general form (in terms of r) 

where ai(l) and bi{l) are constants depending on I. As we can see, these functions have a number of poles 
determined by the zeros of the polynomial in the denominator. If any of these poles happen to lie close to 
our path of integration we can expect large variations in g. This is indeed what happens for I = 7, causing 
convergence problems for the spectral code (which relies crucially on the smoothness of the solution). We 
have tried to use other integration routines and found it straightforward to improve the accuracy of the 
obtained frequencies. However, the main goal here is not to obtain accuracy in this highly unphysical toy 
problem, but rather to assess the pros and cons of our proposed characteristic scheme. Hence, the convergence 
problems displayed in table [T] only serve as a nice reminder of the kind of problems that can be encountered 
in real situations. We should note here that if a pole should happen to lie exactly in the path of integration 
the method presented here will fail. This seems highly unlikely to happen however and is not related to the 
characteristic approach, but rather to the analytic continuation and introduction of a complex r. 

3 Quasi- normal modes of uniform density "stars". 

We now turn to the slightly more relevant problem of finding the axial ly- modes of a constant density star. 
This problem has been treated many times before {e.g. [111 [HI [IHl HI] ) allowing us to compare our results 
to the existing ones. Since the behaviour near null infinity is most strongly dependent on the frequency and 
not on the particulars of the wave function(s) near or in the star we believe that this comparatively simple 
problem serves as a sufficient test of our method as far as non-rotating stars are concerned. In order to 
span as large a region of frequency space as possible we concentrate on an ultra compact model for which 
R — 2.26M, i.e. quite close to the Buchdahl limit R = 2.25M. The w-mode spectrum can then be 
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loosely divided into three parts (see [3]). The trapped modes jl9| are characterised by a low damping rate, 
i.e. a small imaginary part of the frequency. These modes exist due to a peak in the potential near r = 3M 
and correspond, in a loose sense, to almost bound states of the Schrodinger-like wave equation. For higher 
"energies" the modes become less and less bound {i.e. more damped), having increasingly large imaginary 
parts of the frequencies. These are the ordinary, or curvature, w-modes. An infinite number of these modes 
is likely to exist. 

There also exist another branch of modes, the interface, or w/j-modes |23| . These are thought to corre- 
spond to scattering of gravitational waves off the "surface" of the star much like hard-sphere scattering of 
sound waves. These rapidly damped modes have large imaginary parts of the frequency. Previous numerical 
surveys have only found a few (typically two or three) such modes. With our method we are able to probe 
larger imaginary parts and discover many more w//-modes, thus raising the question whether these modes 
also form an infinite family. 

We now proceed to show that the analysis carried out in the preceding section generalises to spherically 
symmetric spacetimes with arbitrary M. We shall only consider the characteristic approach here. It is 
straightforward to show that the wave function has the asymptotic behaviour 



1 - 



il{l + l) 12iMuj + 2l{l + l) + P{l + lY 



2ujr 



^2iuir4iMui 



8uj'r' 
i[16APuj^ - l{l + 1)] 
2ujr 



Using the relation ([7]) we can write 



C{r - 2M) 



(JiiMu 



1 



8iM^Lu 



(32) 
(33) 

(34) 



The outgoing solution approaches a constant whereas the ingoing solution, which is proportional to e^*'^''* , 
decays exponentially for damped modes. The phase function has the leading order behaviour 



il{l+l) ^ 

'),,„2 '-'out 



2iuje 



so that, for damped modes. 



9 ~ 5out 



il{l 



1) 



2wr2 



(35) 



(36) 



making manifest the problem of separating a general solution from a purely outgoing one. Introducing the 
complex r coordinate given in equation (j24p we turn the damping exponentials into simple trigonometric 
functions so that for large g and Ci„ ^ 



For Ch 



2icje2*l'^leCin 

0, i.e. the outgoing solution, we have instead 

il{l + l)uj 
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(37) 



(38) 



It is clear that if the boundary condition 5 = is imposed "exactly at infinity" , there is no contamination 
of the unwanted solution whatsoever. However, one may ask whether this preferred state of affairs prevails 
if one is forced to impose the boundary conditions at some large but finite radius g. Putting the expansions 
equal we find that the asymptotic condition for an outgoing solution to coincide with a generic solution at 
some given radial coordinate g = goo (say) is 



„2i\ 



l{l + l) 



Co 



(39) 
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We see that for any finite there exist complex C\n and Cout such that numerical confusion between the 
wanted Cin = solution and a mixed in- and outgoing solution aris^l- However, it is also clear that this 
solution has Cin ~ (1^^ | Poo Cout and is therefore very close to the wanted solution. Thus, the only situation 
where this can cause worries is when the ingoing solution is rapidly decaying as a function of radius in some 
region (so that when integrating from large radii this solution blows up) . Such a blow up will indeed happen 
in black hole spacetimes and one should also be careful if stars more compact than R < 3M are considered 
(due to the peak in the effective potential). 

3.1 Numerical implementation 

Performing the compactification as described in the previous section, i.e. changing the radial coordinate to 

^^_r__R^r2 ^ r^R + r2^—^ (40) 
r — R + r2 1 + X 

and using the phase function ()22|) we obtain the equation 



2^2 2 , -2. 4r2 fM . \ 2r2 f 6M + 
19 +e z:^^ ( — + e , 1^2 ( 73 ^] (41) 



(.T + l)2^ {x + iy\r^ J'' {x + l)^\ 



for the exterior vacuum perturbations. This equation is solved by the same routine that was applied to the 
toy problem of the previous section, supplying a function gf{oj) = g{x = l,uj). The interior is governed 
by the wave equation (jlSp with the potential given by Since our main aim here is to demonstrate 

the usefulness of the approach to the exterior perturbations we used a very simple code for the interior 
problem. The equations were taken from [24j and are listed in the appendix. We integrated the wave 
function in Schwarzschild coordinates using an off the shelf Runge-Kutta routine. At the surface of the star 
the phase function was evaluated, remembering eq. (jlOp thus giving rise to the ftmction gl{uj) = g{r = R,lo). 
Now gf{uj) and 5^ (w) are guaranteed for every lu to satisfy the boundary conditions at infinity and r = 
repectively. It follows that the quasi-normal mode frequencies are determined by the continuity of g at the 
surface of the star, and are hence given by the roots of 

5r =.9fH-.9fM=0 (42) 

We used the same strategy to find the modes as in the toy problem. In figure [3] we display contour plots of 
\gl°^\. The modes appear as "islands" in this plot and using the locations of these islands as initial guesses 
for our Mtiller root finder we obtain the roots (i.e. the quasi-normal mode frequencies). A sample of the 
obtained frequencies are given in table [31 They are shown in figure [3] as black squares. All results agree 
well with those of Kokkotas [13 US] and exactly with those of Tominaga, Saijo and Maeda [H]. In the left 
panel of figure [3| we show the interface modes. We have managed to locate eight such modes. This raises 
the question if these modes constitute an infinite family. For the most rapidly damped modes in this class 
we encounter convergence problems both in the interior and the exterior routines. The exterior problem is, 
however, less severe than the interior one and a relative accuracy of about one part in 10^ in the eigenfunction 
is achieved with moderate resolution. The interior solution is more difficult to determine accurately. A quick 
glance at the eigenfunctions reveal why, they are rapidly oscillating and exponentially growing suggesting 
that a phase function approach would be beneficial also for the interior problem. As stressed above, we are 
mainly concerned with the methods here and did not consider it relevant to try to improve the accuracy of 
the interior solutions for this paper. From a theoretical point of view it would, however, be interesting to 
settle the question on the number of interface modes. In the right panel of figure [3| we display the trapped 
and curvature modes. For these modes we do not have any convergence problems. 



^By this we do not mean that there is necessarily an exact coincidence of a mixed in/outgoing solution and the purely 
outgoing solution, only that two such solutions are numerically indistinguishable. 
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Table 3: Dimensionless frequencies ui = iVy/BF/3M of the axial I = 2 modes in a i? = 2.26M constant 
density star. The trapped modes, defined to be those whose frequency satisfy SR(a;)^ < 14iax, are marked 
with a Here Knax ~ 0.1513/M^ is defined to be the value of the effective potential in the Regge- Wheeler 
equation at the peak. The interface modes are marked by All digits shown should be correct and in many 
cases the accuracy is much better than displayed. 



11 




12 3 4 



Re((o) Re((B) 

Figure 3: We plot the contours of |(7*°*(d))|, where ui — usy/ /iM. Quasi-normal modes correspond to zeros 
of this function and appear as "islands" in the plot. In the left panel we show the region in the complex Cj- 
plane containing the interface modes. The shaded areas have not been covered due to convergence problems, 
mainly in the interior code, see main text for a discussion. Modes are displayed as black squares. The right 
panel shows the region containing the trapped modes and (a subset of) the curvature modes. 

4 Conclusions 

We have discussed a new approach to the quasinormal-mode problem in general relativity. By combining a 
characteristic formulation of the perturbation equations with the integration of a suitable phase-function for 
a complex valued radial coordinate, we have reformulated the standard outgoing-wave boundary condition 
as a zero Dirichlct condition. This brings a number of important advantages over previous strategies. 
The characteristic formulation permits coordinate compactification, which means that we can impose the 
boundary condition at future null infinity. The phase function avoids oscillatory behaviour in the solution, 
while the use of a complex radial variable allows a clean distinction between out- and ingoing waves. We have 
demonstrated that the method is straightforward to implement, and our analysis of two simple toy problems 
shows that it can lead to high precision numerical results. It is worth noting that the generalisation to 
unstable modes is straightforward, only requiring an alteration of the integration contour. 

Even though the numerical results we have presented are in themselves of no great interest, the new 
method may represent a breakthrough in this area. Most importantly, we have every reason to believe that 
it should generalise to the problem of rapidly rotating relativistic stars. In that case one would no longer deal 
with a simple one-dimensional wave equation in the exterior vacuum. The lack of workable implementations 
of the required outgoing- wave conditions in that problem has been holding back progress for many years. 
Our new approach may be the key that unlocks this problem. Of course, we are still some steps away from 
implementing our new ideas for rotating star spacetimes. Most importantly, we need to consider how the 
different ingredients in our prescription generalise to this more complicated problem. This is an interesting 
challenge and we would hope to make progress on it in the near future. 
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A Constant density sphere — w-modes. 

We test the code by computing the w-modes of a uniform sphere of mass M ^ 0. In the interior the 
background line element can be written 



e^'Wt^ + e^^dr^ + r^{de^ + sin^ 



where 



I 2M 
e- = - I 3a/i - — 



2Mr2 



1 - 



see e.g. 



The perturbation equations can be put in the form 
dr 



dr 



where 



Xi = icjr^^'+^V 
X2 = ^e'''--''uo'^r-^ip 



(43) 

(44) 
(45) 

(46) 
(47) 



(48) 
(49) 



and Vc is the central value of v. The function Lp can in this context be viewed as auxiliary and defined by 
equation 

At the centre of the star the regular solution behave as 



Xi=Xi{l 



e-^-^uj^ + {l + 2)§,{e-''^-2l) 



2(2/ + 3) 



X2=Xi |-(/ + 2)^ 
where Xi is an arbitrary constant 



{I + 4)e2"-co2 - (/ + 2)(/ - l)|^(e-"- + 2Z + 6) ^ 
2(2? + 3) ^ 



- 0{r' 



(50) 
(51) 
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